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Analysis of an Arctic sea ice loss model in the limit of a discontinnons albedo 


Kaitlin Hil0, Dorian S. AbbolU, and Mary SilbeiH 

Abstract: As Arctic sea ice extent decreases with increasing greenhouse gases, there is a growing interest in whether 
there could be a bifurcation associated with its loss, and whether there is signihcant hysteresis associated with 
that bifurcation. A challenge in answering this question is that the bifurcation behavior of certain Arctic energy 
balance models have been shown to be sensitive to how ice-albedo feedback is parameterized. We analyze an Arctic 
energy balance model in the limit as a smoothing parameter associated with ice-albedo feedback tends to zero, which 
introduces a discontinuity boundary to the dynamical systems model. Our analysis provides a case study where 
we use the system in this limit to guide the investigation of bifurcation behavior of the original albedo-smoothed 
system. In this case study, we demonstrate that certain qualitative bifurcation behaviors of the albedo-smoothed 
system can have counterparts in the limit with no albedo smoothing. We use this perspective to systematically 
explore the parameter space of the model. For example, we uncover parameter sets for which the largest transition, 
with increasing greenhouse gases, is from a perennially ice-covered Arctic to a seasonally ice-free state, an unusual 
bifurcation scenario that persists even when albedo-smoothing is reintroduced. This analysis provides an alternative 
perspective on how parameters of the model affect bifurcation behavior. We expect our approach, which exploits the 
width of repelling sliding intervals for understanding the hysteresis loops, would carry over to other positive feedback 
systems with a similar natural piecewise-smooth limit, and when the feedback strength is likewise modulated with 
seasons or other periodic forcing. 
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1 Introduction 

A current area of interest in climate science is how perturbations, such as increased greenhouse 
gases, will affect Arctic sea ice, and what will be the manifestation of this change. One possibility 
is that the ice, which is currently present year-round, and which cycles between different extremes 
with the seasons, will cross a “tipping point” to a completely different state. Mathematically, a 
tipping point may be associated with a bifurcation, after which further parameter changes result 
in a large and sudden change in system behavior that is not easily reversed (i.e., a hysteresis loop 
forms). In light of the rapid decrease in summer sea ice extent in the Arctic region during the past 
decade m, a potential “tipping point” in summer sea ice has been considered in recent studies 
of conceptual dynamical systems models [H El 0 Ell. Tipping points have also been extensively 
investigated in large-scale global climate models, which examine the loss of sea ice from the Arctic 
region in the coming decades under various greenhouse gas emission scenarios PEHESlEe]. 

In this paper we analyze a single-column energy balance model for the Arctic region, originally 
proposed by Eisenman and Wettlaufer [7] (hereafter EW09) as a simple conceptual framework for 
investigating the role sea ice thermodynamics may play in bifurcations associated with the loss of 
Arctic sea ice. Their model took the form of a nonautonomous, one-degree-of-freedom ordinary 
differential equation, which captured the yearly variation in solar forcing, averaged over the Arctic 
region, as well as seasonal variations in heat transport from lower latitudes into the Arctic. It 
incorporated the classic “ice-albedo feedback,” a positive feedback mechanism that can lead to 
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coexisting stable states: one state represents a perennially ice-covered Arctic with high albedo 
(i.e., high reflectivity), and the other a perennially ice-free Arctic with low albedo [2T]. Their 
model also incorporated the sea ice self-insulation feedback, a thermodynamic feedback mechanism 
where thicker ice insulates itself better. EW09 used their model to show how sea ice self-insulation, 
which allows thin ice to grow more rapidly than thick ice, may mitigate, to a degree, positive ice- 
albedo feedback effects. Their model shows a smooth transition, under increases in atmospheric 
greenhouse gases, from a stable perennially ice-covered Arctic to a stable seasonally ice-free state. 

The model in EW09 evolves an energy density E{t) that is measured relative to an ocean mixed 
layer at the freezing point; when E < 0, E represents the average Arctic sea ice thickness, and when 
E > 0 it represents a mean temperature of the Arctic ocean mixed layer. As the energy fluxes 
change with the state of the system (ice vs. ocean), this conceptual framework leads naturally 
to a formulation of the problem as a piecewise-smooth (PWS) continuous dynamical system |24j : 
across the “discontinuity boundary” at E = 0, the vector held is Lipschitz continuous but not 
differentiable. However, in the limit where we take the albedo function to be piecewise-constant, 
the vector field is generally not continuous across this boundary. Incorporating a piecewise-constant 
albedo function is a standard approximation exploited in analytic studies of energy balance models 
first developed to investigate planetary energy balance; in such studies the albedo jumps between 
high and low values at the latitude of the ice cap boundary [a [201 [23]. In our investigation a 
boundary in phase space is crossed when the sea ice is lost and then again when it is recovered over 
the course of the yearly cycle of a seasonally ice-free solution. We focus our analysis on the key 
role this discontinuity boundary plays in the bifurcation structure of the problem. For a review of 
PWS systems and bifurcations involving discontinuity boundaries, see [6|. 

EW09 modeled the transition between when the Arctic is ice-covered and ice-free by introducing 
a smoothing energy scale parameter AE over which the albedo transitions smoothly between a high 
value (for sea ice) and a low value (for open ocean). The smooth albedo transition is meant to 
parameterize the spatial variation of ice and capture the albedo dependence on ice thickness [7]. 
However, model results are sensitive to the size of this energy range: transitions between solution 
states, via pairs of saddle-node bifurcations, are easily “smoothed out” by increasing the energy 
range AE in the albedo transition [Tj. We take this parameter to zero to formulate our system. An 
immediate consequence of removing AE from the albedo function is that the resulting system then 
admits “sliding regions,” which classify it as a Filippov system [9|. In general, a Filippov system 
is one for which the flow components normal to the discontinuity boundary at a point may have 
opposite signs on either side of the boundary. When this occurs a trajectory may “slide” along the 
boundary through so-called sliding regions [6|. Specifically, if we introduce time as a dynamical 
variable to render the model of EW09 an autonomous system, then it takes the form 

('G+(r,E), 
dt \G_(r,E), 


Here G+ and G- are period-one functions of r, which is measured in units of years. E = 0 
represents a discontinuity boundary in the (r, E)-phase plane. Note that the dynamics in the 
E = 0 boundary have not been fully specified in the above formulation; this ambiguity in the 
dynamics is inconsequential for evolving a solution through the discontinuity boundary provided 
that dE/dt has the same sign as E —0^ as for E —>■ 0“. However, for the model of EW09 in the 
Filippov limit, there are necessarily intervals of r for which G+(t, E) and G_(t, E) have opposite 
signs as E —>■ 0+ and E —>■ 0“, respectively. These r intervals are associated with “repelling sliding 


E > 0, 
E < 0, 


( 1 . 1 ) 

( 1 . 2 ) 
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intervals” in the model, as dE/dt > 0 for ill —)• O'*' and dE/dt < 0 for iH —>■ 0“. Uniqueness of 
solutions with initial conditions in such sliding intervals is manifestly violated [6]. Solutions that 
enter repelling sliding intervals have the property that they may slide, for an indeterminate time, 
before either jumping off the discontinuity boundary in either direction or reaching the end of the 
sliding interval. (Note that the maximum amount of time a solution may slide is equal to the width 
of the sliding interval, since solutions E{t) evolve r according to dr/dt = 1.) In this way the lack 
of uniqueness leads to families of solutions that interact with the sliding intervals. These families 
of nonunique solutions, when omitted from the analysis, create “gaps” in the bifurcation diagram 
where solution branches become discontinuous. 

Our analysis contributes a new case study of how the Filippov limit of a dynamical system, 
where the vector field is discontinuous and which here is easier to analyze in closed form, can be 
utilized to determine properties of the original continuous system. The qualitative comparison 
between a smooth system and its nonsmooth counterpart has been of significant recent interest in 
the context of applications uniiiT] as well as in general mathematical investigations that consider 
either “smoothing out” a PWS system [I2lll3l[25] or transforming a smooth system into a PWS 
one [5]. We explore how the relative widths of the sliding intervals affect the bifurcation structure, 
which provides insight into how model parameters impact the transition from an ice-covered to 
an ice-free Arctic. For example, we show that a small hysteresis loop, involving seasonally ice- 
free states and associated with the loss of perennial ice, can be moved along the seasonal solution 
branch by varying the relative widths of the sliding intervals. In particular, it can be moved so 
that it is associated with the loss of the perennially ice-free Arctic state rather than the perennially 
ice-covered Arctic state. Additionally, we use insight gained from this perspective to determine 
regions of parameter space where the transition from an ice-covered to a seasonally ice-free Arctic 
is large. 

This paper is organized as follows. Section [2] describes the original formulation of the Eisenman 
and Wettlaufer model [7], as well as variations on it BE], and introduces our general Filippov 
formulation of the model. In Section [3] we derive, in turn, the existence and stability conditions for 
the three types of periodic solutions: perennially ice-free, perennially ice-covered, and seasonally 
ice-free. Section 0] presents a case study using the model given in [8|. In this case study we 
compare numerical bifurcation diagrams for the albedo-smoothed system with analytical bifurcation 
diagrams for the corresponding Filippov system, and we explore the various parameters and how 
they impact the bifurcation structure in an alternative manner to previous parameter studies of 
the model. A summary and discussion of our results concludes the paper in Section [5l 


2 Model Description 

The model developed in EW09 is an energy balance model, where the energy per unit surface area, 
E, is measured relative to an Arctic ocean mixed layer at the freezing point. An increase in the 
energy density relative to this baseline leads to warming of the mixed layer to a temperature T (x E 
above the freezing point, while a decrease in the energy density leads to sea ice growth to an ice 
thickness hi oc —E. Specifically, the state variable E in EW09 has the physical interpretation that 

f -Lihi, E <0, , ^ 

[CsT, E>0, 

where the proportionality constants, Li and Cg, correspond to the latent heat of fusion for sea ice 
and the ocean heat capacity per unit surface area, respectively. 
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The evolution equation for E{t) in EW09 is derived by considering heat flux balances at the 
surface, and it takes the form [7] 


— = (1 - a{E)) Es{t) - {Eo{t) + FT{t)T{t, E)) + Eb - uo'H{-E)E. (2.2) 

at 


The energy source term (1 — a{E))Fs{t) represents the absorbed incoming radiation, where Fs{t) 
is the solar radiation flux, averaged over the Arctic region, and (1 — o.{E)) is the fraction of the 
incoming radiation absorbed {a{E) is the albedo). The energy loss term (Tb(t) + FT{t)T{t, E)) 
represents outgoing radiation, which increases with the surface temperature T{t,E)-, it has been 
linearized about the freezing point T = 0. The flux terms Fs{t), Tb(t), and FT{t) are all periodic 
with period 1 year. Here Fb represents a small, constant basal heat flux into the Arctic mixed 
layer/sea ice from the deep ocean. The term i/qH{—E)E captures sea ice transport out of the 
Arctic, where 7^(x) is the Heaviside function; this term is only present if there is sea ice. In EW09, 
the albedo function smoothly transitions from its ice value (ctj) to its ocean mixed layer value (ami) 
and is given by 

/ 771\ Olml + ai ami ~ 1. ^ \ /'o 

o{S) = tanh j . (2.3) 

where AE is the previously mentioned energy scale for smoothing. 

Sea ice thermodynamics enters the model in EW09 by the surface temperature function T(t, E), 
which depends on the state of the ice/mixed layer: 


T{t,E) 


'eICs, 

< 0 , 

/ (l-«i)Fg(f)-Fo(t) \ ( E \ 

A FT(t) ) \E-kiLi/FT{t) ) ' 


E > 0, 

E < 0, (1 - ai)Fs(t) > Fo{t), 
E < 0, (1 - ai)Fs(t) < Fo{t), 


where ki is the thermal conductivity of sea ice. Note that when E > 0, no sea ice is present, 
so T is determined by (Et]). When E < 0, sea ice is present, and T depends on whether or not 
the ice is melting at the surface. When the outgoing radiation exceeds the incoming radiation 
((1 — ai)Fs(t) < Eo(t)), surface temperature depends on the thickness of the sea ice, hi oc —E, 
and it is determined, in a large Stefan number approximation, from surface flux balance, fcjT /hi = 
(1 - ai)Fs{t) - Fo{t) - FT(t)T. (See [7] for more details.) On the other hand, when E < 0 and 
incoming exceeds outgoing radiation, the ice is melting, and in a large Stefan number approximation 
the surface temperature immediately warms to the melting point, defined here to be T = 0. 

A variation of the model in EW09, due to Abbot, Silber, and Pierrehumbert [T], replaces the 
seasonally varying Eo(t) with a state-dependent model and incorporates possible effects due to cloud 
cover into the albedo function. In another variation, Eisenman [8] simplified the EW09 model by 
introducing a sinusoidal approximation for each of the seasonally varying terms Fs(t) and Eo(t), 
while neglecting any seasonal variations of FT{t) in (12.2p . This is the version of the model that we 
use in our case study in Section 01 

Our analysis is based on these Arctic sea ice models in the following form: 


^=E(r,E)-Er(r,E), 

- = 1, 

dt ’ 


(2.4) 
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where after suitable nondimensionalization of [ 8 ], the surface temperature is given by 


T{t,E) 


E, 


< 0 , 

F 

B 

\ 



E>Q, 

E <0, F >0, 
E <0, F <0, 


(2.5) 


where C > 0 is a dimensionless thermodynamic parameter associated with sea ice. This is the form 
of (| 2 . 2 p from EW09 in the case that basal heat transport and sea ice transport out of the Arctic are 
neglected, i.e., Fb = f = 0. Moreover, we have neglected any temporal variation of FT{t) in (12.211 by 
replacing it with the constant B. (Similar approximations were made in some numerical simulation 
results in [S].) Note that although T{t,E) is piecewise-defined, it is Lipschitz continuous, so it does 
not contribute to the discontinuity in the vector field at E = 0. 

In the limit of no albedo-smoothing, F (r, E) in (j2.4p is given by 


- Ei^+{t) = F+{t), E>0, 
^ \ {l-Aa)Es{T)-Ei,.{T) = F_{T), E<0, 


( 2 . 6 ) 


where F±(t) are (nondimensionalized) period-one functions of r. The albedo function here is the 
limit as AE ^ 0 of the albedo function after nondimensionalization in [ 8 ]; the term 


Aa = 


2 (1 — (aj -|- ami)l‘^) 


e[o,i] 


(2.7) 


represents the difference between the ice and ocean mixed layer albedos, scaled as part of the 
nondimensionalization process by twice the mean coalbedo (coalbedo = 1 - albedo), which is 
approximately 1. The fluxes Es{t) and Ei^±{t) characterize incoming shortwave and outgoing 
longwave radiation terms, respectively, akin to Fs{t) and To(t) in (|2.2p . The subscript ± on the 
outgoing longwave radiation term indicates that it may depend on whether there is ice in the 
Arctic or not (e.g., as in the variations of the model of EW09 investigated by Abbot, Silber, and 
Pierrehumbert [T]). In our case study in Section [5] we use 


E) = < cos(27rT)) - {Lm + La cos(27r(T - (/>))), E > 0, 

|E_(r) = (1 - Aq;)( 1 - 5'aCOs(27rT)) - (Lm + TaCOs(27r(T - (/>))), E < 0, 

which was generated by taking the limit as AE —>■ 0 of E(r, E) from [8], 

E(t, E) = (1 -|- Aa tanh(E/AE))(l — Sa cos(27rr)) — {Lm + La cos(27r(r — 0))). (2.9) 


This represents a sinusoidal approximation to the periodic functions Fs{t) and Eb(t) from EW09 
[ 8 ]. We note that Lm, the annual mean outgoing longwave radiation when T = 0, decreases as 
greenhouse gas concentrations increase [ 8 ] and is the parameter we vary in our bifurcation studies. 
The general system we consider can be summarized by combining (I2.4I) - (I2.6I) : 


dE 

dt 


'E+(T)-EE(r), 

< E_(r), 

SF.{t)/{C-E{t)), 



E > 0, 

E < 0, E_ > 0, 
E < 0, E_ < 0, 


( 2 . 10 ) 


We will refer to this system as the “Filippov system,” and we will refer to the original system with 
albedo smoothing, i.e., AE 7 ^ 0 in (12.9h . as the “albedo-smoothed system.” There are three possible 
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(a) 



ta ib 


(b) 



T 

(c) 


Figure 1: (a) Sample plot of F+(r) (solid red curve), which applies for E > 0, and F-{t) (dashed 
blue curve) for E < 0, from (|2.6p over one period of the seasonal forcing in EW09 [7], after 
nondimensionalization. Repelling sliding intervals are shaded; triangles indicate the point where a 
trajectory may enter a sliding interval. For perennially ice-free solutions, ta is when the minimum of 
E{t) occurs as E —> O"*". For perennially ice-covered solutions, ice ablation initiates and terminates 
at tb and tc, respectively. For seasonally ice-free solutions, we show in Section [3] that ta is the last 
possible freezing time and tc is the last possible melting time. (b),(c) Blow-up of the direction field 
for the nondimensionalized EW09 model near the sliding intervals; nullclines are given in red for 
E > 0 and blue for E < 0. Shaded regions represent the sliding intervals = [ta,tb], associated 
with “spring,” and S 2 = [tc,td\, associated with “autumn.” 


discontinuity boundaries in this system: one at E = 0 that is a result of the albedo jump, and two 
corresponding to roots of E_(t), which are associated with changes in the thermodynamics between 
growing and melting sea ice. Across each E_ = 0 boundary the slope of solution trajectories jumps 
for E < 0 solutions. For our purposes, though, these boundaries associated with F- = 0 are 
benign, in the sense that trajectories cross them with transverse speed given by dr/dt = 1. In 
contrast, the discontinuity boundary E = 0 associated with the albedo jump may introduce sliding 
regions; our analysis will focus on how solutions interact with these. Solution behavior on either 
side of the discontinuity boundary at E = 0 is characterized through F±{t), given by (12.6p . The 
time intervals for which E+(t) > 0 and E_(r) < 0 correspond to repelling sliding intervals for 
the Filippov dynamical system (|2.10p . These time intervals, one in spring and one in autumn, are 
indicated by shading in Figure [H Note that the thermodynamic parameter ( and the parameter 
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B from (12.101) do not enter into determining the intervals of sliding. The discontinuity boundaries 
meet at values where E = 0 and F_ = 0 (e.g., the points (4,0) and (4,0) in Figure [I]). Since the 
F- = 0 boundaries are “crossing regions” they do not interfere with the crossing or sliding behavior 
on the E = 0 boundary at these points. 

Note that the only way for a solution trajectory to enter a repelling sliding interval is if i? = 0 
at the start of a sliding interval (e.g., if FI —>• O"'' as r —>■ ta or if Fi —> 0“ as r —> tc in Figured]). A 
technique called “Filippov’s convex method” can be used to define the behavior of the solution in 
the boundary when it is sliding [9]. This method can in general reveal diverse phenomena such as 
equilibria of the sliding flow [6], but in our case it produces trivial sliding motion {dE/dt, dr/dt) = 
(0,1) due to the system dynamics dr/dt = 1. 

We do not consider the possibility of “attracting sliding intervals” for this model since the 
primary source of discontinuity in the vector field is due to the albedo jump Aa > 0 at Fi = 0, 
i.e., through the contribution (l±Aa)F'<j(T) > 0 to (|2.6)1 . This contribution is enough to ensure that 
F+ir) > F- (r) over most of the yearly cycle of forcing. An attracting sliding interval would require 
that dE/dt < 0 for Fi —O’*" and dE/dt > Q for E ^ , which would require F4(r) < Fh(r) over 

some r-interval. Attracting sliding intervals do occur, for some parameter sets, in the Filippov limit 
of the simplihed model investigated in [8] , but they are a nonphysical consequence of the sinusoidal 
approximations used for the periodic forcing functions. For this reason we focus our bifurcation 
analysis in Section 0] on parameter sets for which attracting sliding does not occur. 

3 Nonsliding Periodic Solutions 

In this section we determine the existence and stability properties of the three distinct types of 
periodic solutions of the Filippov system ()2.10p : 

A. Perennially ice-free {E{t) > 0 for all r). 

B. Perennially ice-covered {E{t) < 0 for all r). 

C. Seasonally ice-free {E{t) takes on both positive and negative values). 

By virtue of their definitions the perennially ice-free (FI > 0) and perennially ice-covered {E < 0) 
solutions avoid the discontinuity boundary E = 0. Nonetheless, the parameter space existence 
boundary for these solutions is determined by the condition that the minimum (maximum) value 
of E limits to Fi = 0. This limiting minimum value of Fi = 0 occurs for perennially ice-free solutions 
at T = ta, while the maximum value of Fi = 0 occurs for perennially ice-covered solutions at r = tc] 
see Figure [I](b),(c). In contrast the seasonally ice-free solutions necessarily cross the discontinuity 
boundary E = 0 twice during their yearly cycle. In Section 13.31 we show that in addition to saddle- 
node bifurcations along the seasonal solution branch, these states may be destroyed, as parameters 
vary, when they collide with the sliding intervals at either (r, Fi) = (ta,0) or {t,E) = (tc,0). 

3.1 Perennially Ice-Free Periodic Solntions 

In this section we derive conditions for the existence, stability, and uniqueness of perennially ice- 
free periodic solutions that do not enter the sliding interval. Specifically, we show that these 
solutions are always unique and stable when they exist, and that this solution branch arises via a 
“grazing-sliding bifurcation” [6| with the discontinuity boundary E = 0 at r = ta- A grazing-sliding 
bifurcation occurs when a trajectory that resides entirely on one side of the discontinuity boundary 
is continuously perturbed into a trajectory that grazes the boundary at a point and, with further 
perturbation, may begin to slide [6|. 

Perennially ice-free solutions satisfy Fi > 0 for all r, so that based on (|2.10p . E evolves according 
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dr 

(For Section [3] we use the notation ^ for simplicity, since ^ = 1.) From this we can derive a 
period-one Poincare return map, which takes an initial condition Eq = F'(to) time tg ^ [Oj 1) to 
its value one period later: 

E{tQ + l) = e ^[E{tQ)+e to + 1)]; (3-2) 

where ^ 

/+(to,r)= [ e^^E+{t)dt. (3.3) 

Jto 

(At this stage we have not chosen a preferred phase to for our map; it is chosen below.) Periodic 
solutions satisfy Fi(to + 1 ) = E(to) = Eq, which determines the unique initial condition for the 
periodic state: 

E*o = (e^ - l)“^-^'°/+(to,to + 1). (3.4) 

The stability of this periodic state is determined by the eigenvalue of the linearized Poincare map, 
which is the Floquet multiplier 

H = dE{to + l)/dE{to) = e~^ . (3.5) 

Since i? is a positive constant, |/r| < 1 and the periodic solution is always stable when it exists. 

Finally, we determine the existence condition for perennially ice-free solutions, which is that 
E{t) > 0 for all values of r E [toj^o + !)• To do this we consider the time at which E achieves its 
minimum, tmin, and let E{tmin) 0 . The conditions for a minimum, that ^ = 0 and > 0 , 
become simply E+{tmin) = 0 and ^ > 0 for E{tmin) 0 , so tmin ia- (See Figure [I](a) for 
a sample plot of T+(t).) Hence we may set to = tmin = ta in (I3.4jl so that perennially ice-free 
solutions exist provided Eq > 0, i.e., provided 


I+{ta,ta + l) > 0, (3.6) 

where ta is determined from 

E+{ta)=0, ^{ta)>0. (3.7) 

dr 

In this analysis of the perennially ice-free periodic solutions we have avoided the sliding intervals 
of the Filippov system by insisting that E{t) > 0 for all r. However, the existence boundary for 
these solutions in parameter space corresponds to conditions where E{ta) —)• 0 , which is also the 
unique entry point to 5i, the springtime sliding interval in Figure [l](b). Thus the bifurcation that 
leads to the creation of the perennially ice-free solutions corresponds to a grazing-sliding bifurcation 
[ 6 ]. In relation to the classification of grazing-sliding bifurcations in [H], this bifurcation is similar 
to the generic case TC 2 , where a sliding periodic trajectory and a nonsliding periodic trajectory 
collide and annihilate each other. 

3.2 Perennially Ice-Covered Periodic Solutions 

For these periodic states, E{t) < 0 for all r. We introduce a time tb E [0,1) when surface ablation 
initiates so that T = 0 in (j2.4p . and another time tc E {tb,tb + 1) when surface temperatures next 
drop below the freezing point. Our focus is on periodic solutions, and we consider those over the 
time interval r E [tb,tb + 1]. There are thus two phases associated with this solution: a summer 




phase, r € [tb,tc\, when ice is ablating from the surface and the surface temperature is fixed at 
T = 0, and a winter phase when T drops below the freezing point, which occurs for r € + !)■ 

Specifically, tb and tc are defined by the conditions 


F_(4) = 0, -^(4)>0, 

dF_ 

F_(tc) = 0, -^(ic)<0, 

as shown in Figure ma). Based on these two phases, ()2.10p simplifies to 


(3.8) 


^ _ f t G [tb,tc], 

dr I > Te{tc,tb + 1), 


(3.9) 


Note that the times tb and tc separate the phases when the perennial ice is melting {E'{t) > 0 
for r G {tb,tc)) from the phases when it is growing {E'{t) < 0 for t G (fc,4 + 1)), since —E is 
proportional to ice thickness when < 0. It follows that the minimum ice thickness over the yearly 
cycle occurs at time t = tc- 

Letting Eb = E{tb) be the initial condition and Ec = E{tc), we find the period-one Poincare 
return map, given implicitly by 

E{h+l)-E^!!^=E,-!^ + {F_), ( 3 . 10 ) 

where 


Ec — Eb + I-{tb, tc), 

{E_)= [ F_(r)dr, 

Jo 

(3.11) 

(3.12) 

I-{tb,tc)^ [ E_{T)dT. 

Jtb 

(3.13) 

The initial condition for a periodic solution, E{tb) = E{tb + 1) ^ E^, is then 

C{E-) I-{tb,tc) 

^ I-{tb,tc) 2 

(3.14) 


We obtain the condition for a perennially ice-covered periodic state to exist by insisting that 
E < 0 for all r G [tb,tb + 1). Since the minimum ice thickness occurs at r = tc, we obtain the 
existence condition as 

E: = (F_) + < 0, (3.15) 

where E* is the value of E{tc) for the periodic solution beginning at E^. We can confirm that this 
solution is stable by computing the Floquet multiplier from the map (I3.10p . We find 


dE(tb 1) 
dEb 


E,=El 


C-E* 

C-Ef 


(3.16) 


Since E^ < E* < 0 and C is positive, the Floquet multiplier must be in the interval (0,1), and 
hence the solution is stable when it exists. 

Again, we have avoided the sliding intervals of the Filippov system by insisting that E^ax = 
E* < 0. However, the existence boundary in parameter space corresponds to where E{tc) 0, 
which is the unique entry point to S 2 , the autumn sliding interval in Figure [11(c). Thus the 
bifurcation that leads to the creation (or annihilation) of the perennially ice-covered solutions also 
corresponds to a grazing-sliding bifurcation. 
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Time 

^ T <C t^n 


tm<T <tf 


tf <T <tb + l 

Equation 

II 

o 

II 

-to 

^ = F^-BE 

o 

II 

-to 

dE (F. 

dr ~ C-E 

State 

T = 0, E <0 

E = T >0 

E,T <0 

Ice Phase 

Ice melting 

No ice 

Ice growing 


Table 1: The three thermodynamic phases of a seasonally ice-free periodic solution, and the 
piecewise-defined differential equation that governs the solution. 


3.3 Seasonally Ice-Free Periodic Solutions 

For this state there are three thermodynamic phases, each governed by a different differential 
equation, that are fully summarized by Table [TJ We take as a starting time for the periodic 
solutions tb, when the ablation phase begins; our initial condition is Eb = E[tb). We also introduce 
a “melting time” tm when the ice-free phase starts, and a “freezing time” tf when the ice-free phase 
ends. The time tb is defined implicitly as in (j3.8j) (see Figure [11(a)), while tm and tf are defined by 


E{tm) = 0, 

(It 


(3.17) 


and are not known a priori. These times are by-products of our analysis, and they determine the 
fraction of the year that this seasonal state is ice-free, tf — tm- Note that our requirements that 
dE/dr > 0 at T = tm and dE/dr < 0 at t = tf ensure that this solution does not enter either 
sliding interval at the E = 0 boundary. This follows from the property that the sliding intervals 
are repelling, so they may only be entered if dE/dr = 0 at their left boundaries (see Figured]). 

Starting with Eb = E{tb) and integrating the governing piecewise-defined differential equation 
through each of the phases indicated by Table dl we find that 


Eitb -\- 1 ) 



(— I- {tbi tm)) ) 


(— I+{tm, tf)) , 


(— ^-{tf, tb + 1)) . 


(3.18) 


From these we determine that a periodic solution Eb = E(tb + 1) ^ E/ <0 exists, provided we can 
find times tm and tf which solve the following implicit equations: 


0 = I+{tm,tf), (3.19) 

0 = I-{tb,tm) -^ I-{tf,tb + !)■ (3.20) 

Moreover, any solution {tm,tf) to this pair of equations ensures existence of a periodic seasonal 
state only if it lies in an appropriate domain determined by the following additional constraints, 
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related to the timeline of Table [TJ 


tb < < < tb +1) 

F_(r )>0 VrG( 4 ,tm], (3.21) 

F+(t/) < 0 , 

where F-{tb) = F^{tb + 1) =0. Other conditions suggested by Tabled! that T_(r) < 0 for 
r G [t/,tb + 1) and F^{tm) > 0, are guaranteed by those included in (I3.2ip (see Figure [11(a)). Once 
the times tm and tf satisfying ()3.19p - ()3.2ip are determined, the initial condition for the periodic 
seasonally ice-free state is then given by the explicit expression in (I3.18p . 

Unlike the perennially ice-free and perennially ice-covered states, the seasonally ice-free periodic 
solutions are not necessarily unique [iiiaii], and they need not be stable. We compute the Floquet 
multiplier dE{tb + l)/dEi,, evaluated at by implicitly differentiating (I3.18P with respect 

to Eb, and find 


dEitjj T 1) 
dEb 


— g ^rn) 

E,=El 


C 




* 

b 


^+i^f) F-{tra) 


Saddle-node bifurcations of the seasonally ice-free periodic solutions occur whenever 


(3.22) 


^ F_{tf) E+{tm) ^ 
\C-E*bJ F^itf) E.{t^) ^ ■ 


(3.23) 


Note that if tf ^ ta or tf ^ td, or if tm —>■ h or tm —>■ tc, then the Floquet multiplier would 
diverge since F^{ta) = F^{td) = 0 and T_(tb) = E-(tc) = 0 (see Figured]). This is relevant, as we 
will see further in this section and in Section near certain existence boundaries for the seasonal 
solutions, and causes them to come into existence as highly unstable solutions when they arise 
through grazing-sliding bifurcations. Additionally, note that the dimensionality of this problem 
precludes the possibility of a negative Floquet multiplier and period-doubling bifurcations, at least 
insofar as any such solutions avoid the sliding intervals, where the uniqueness of solutions breaks 
down. 

Finally, we may determine the seasonal state existence region in parameter space by considering 
all possible mechanisms for creation of seasonal ice-free solutions: 

1. Saddle-node bifurcations for which (I3.23p is met. 

2 . Solutions (tm,tf) to (|3.19p . (|3.20l) entering through the boundaries of (tm, tj)-space associated 
with the inequalities (I3.2ip . which are tm = tfe, tm = to tf = td, and tf = ta + f- These 
follow from the conditions E-{tm) > 0 and E^{tf) < 0 in (|3.21l) . Below we argue that {tm,tf) 
solutions to (I3.19h . ()3.2np may in fact only enter or leave the {tm,tf) region of existence at 
the boundaries tm = tc and tf = ta + f, which correspond to solutions that enter a sliding 
interval and hence are associated with grazing-sliding bifurcations. 

To eliminate the possibility of (tm,t/) solutions to (I3.19l) . (|3.20p emerging from the boundaries 
tm = tb and tf = td, we consider the relationship between these times and the sliding intervals of the 
discontinuity boundary, as illustrated by Figured) The only way a solution can have a zero at either 
tb or td is if it enters the sliding interval at ta or tc, respectively. However, since the vector field 
at T = ta points downward for E —)• 0“(see Figure [11(b)), tb cannot correspond to a melting time. 
By a similar argument, td cannot be a freezing time (see Figure dKc)). Hence {tm,tf) solutions do 
not emerge from the boundaries tm = tb and tf = td- The possibility of {tm,tf) solutions emerging 
from the boundaries tm = tc and tf = ta +1 is confirmed through the considerations above, namely, 
that ta could correspond to a freezing time and tc could correspond to a melting time. Note that 
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these limiting values are associated with the onset of the springtime and autumn sliding intervals, 
respectively. Thus, ta can be considered the last possible freezing time and tc the last possible 
melting time. 

To summarize, in addition to saddle-node bifurcations, seasonal solutions may come into exis¬ 
tence in parameter space when we can find which satisfy 

0 = I+{tc,tf), 

0 = tc) H-^—I- I-{tf, h + 1) 

at the tm = tc boundary, or 


0 — I+{tm, ta + 1 ), 

0 = /_ (4, + I- {ta, h) 


(3.25) 


at the boundary where tf = ta + 1- These solutions come into existence via grazing-sliding bifur¬ 
cations and are similar to the generic case TCi in [Hj, although with repelling sliding instead of 
attracting sliding. 


4 Case Study Bifurcation Analysis 

In this section we develop a case study comparing bifurcation results from the nondimensionalized 
albedo-smoothed model given by Eisenman [8] with those associated with its Filippov limit achieved 
by taking the parameter AE —>■ 0. We demonstrate the similarity of analytical results in the 
Filippov limit to numerical results of the corresponding smoothed system from [8] through selected 
bifurcation diagrams. The section concludes with an exploration of how varying the relative widths 
of the two sliding intervals in the Filippov system affects the bifurcation diagram. 


4.1 Characterization of Bifurcations 

The simplified Eisenman model is given by (12.4p . (j2.5jl . with sinusoidal approximations to the 
periodic forcing functions [8]: 

F(r, F) = (1 -|- Aatanh(F/AF))(l — Sa cos(27rr)) — {L^ + La cos(27r(r — </>))). (4-1) 


Here, Aa > 0 is proportional to the jump in albedo between sea ice {E < 0) and open ocean 
{E > 0) as defined in (j2.7jl . AE is the energy scale for smoothing, Sa > 0 captures the seasonal 
variation of incoming shortwave radiation. La > 0 measures the seasonal variation of longwave 
radiation, cj) E controls the phase shift between the peaks in incoming shortwave and 

outgoing longwave radiation, and Lm > 0 is the mean of the outgoing longwave radiation when 
T = 0. Note that in [8], (/> E [0,1); we have chosen to center the interval at </> = 0. In the Filippov 
limit F{t, E) takes the form 


F{t,E) 


F+(r) = (1 Aq;)( 1 - S'a cos(27rT)) - {L^ + La cos(27r(T - (/>))), E > 0, 
E-(t) = (1 - Aq;)( 1 - 5'aCOs(27rT)) - (L^ + La cos (27r(T - 0))), F < 0. 


Figure [2] shows examples of the three types of periodic solutions constructed from the analytical 
results of Section [3l using the case study model (|2.101) . (14.21) for parameter sets where the solutions 
exist stably. 
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Perennially ice-free 


Seasonally ice-free 
Perennially ice-covered 


Figure 2: Examples of stable periodic solutions of (l2.1Up . where F±(t) are given by (j4.2ll . for the 
default parameters Sa = 1-5, La = 0.73, (f> = 0.15, B = 0.45, ( = 0.12, and Aa = 0.43, as in Table 
1 of [8]. We set = 0.92 for the perennially and seasonally ice-free trajectories and = 1.1 for 
the perennially ice-covered trajectory. 


In our bifurcation analysis, as in previous studies PEIIISI [H] of the model introduced by 
Eisenman and Wettlaufer [7], we consider a decrease in the mean value of the outgoing longwave 
radiation (for T = 0) to serve as a suitable proxy for increases in greenhouse gas emissions. For 
example, the mean value of Fb(0 ™ (12.21) . which is denoted Lm after nondimensionalization in 
(j4.2h . is expected to decrease approximately logarithmically with increases in atmospheric CO 2 
concentration [22]. Thus, increases in greenhouse gases have the effect of decreasing the widths of 
the sliding intervals in Figured) The default present-day value of is given as 1.25 in [8| and is 
estimated based on observations mM- 

The Eisenman model (|2.4h . (|2.5h . (|4.1h captures the same bifurcation behavior as the model in 
EW09, as well as the results of several other single-column energy balance and global climate models 
[8]. Through extensive numerical simulations for various parameter sets, Eisenman identihed four 
bifurcation scenarios, representing the range of results from these previous models, which also occur 
in his model. Eor the purposes of our analysis, this model has the benefit that in the Filippov limit 
the integrals associated with solutions, derived in Section |3l can be evaluated in closed form. On 
the other hand, for some parameter sets the model violates the assumption that F+(t) > F_(r), 
for all T, made throughout previous sections. In this analysis we will not consider parameter sets 
for which this violation leads to attracting sliding intervals. 

An example bifurcation diagram of the Filippov system is shown in Figure[3j(a). Two bifurcation 
diagrams of the albedo-smoothed system are given for comparison, in (b) with a small amount of 
smoothing (AE = 0.02), and in (c) with the “default” amount of smoothing used in the analysis 
of [8| {AE = 0.08). We make the following remarks about Figure [3) 

1. Since decreases with increased greenhouse gases, we have chosen to show decreasing 
from left to right in the bifurcation diagram (so that greenhouse gases increase from left to 
right). Solutions progress from the current state of a perennially ice-covered Arctic, through 
a seasonally ice-covered state, to a state in which the Arctic is perennially ice-free at high 
greenhouse gas levels. The transition from seasonal ice to a perennially ice-free state results 
from a saddle-node bifurcation on the seasonal branch. 

2. Notice that we have chosen to represent the periodic solution branches in the bifurcation 
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(a) 




(b) (c) 


— Perennially ice-free 

— Seasonally ice-free 

— Perennially 
ice-covered 


Figure 3; Bifurcation diagrams for (12.41) . (12.511 . where F{t,E) is given by (14.ip . with (a) AE = 0, 
(b) AE = 0.02 (colored curves) and AE = 0 (gray curves), and (c) AE = 0.08 (colored curves) 
and AE = 0 (gray curves), for the default parameters, as in Figure [2l In (a) the diagram on the 
right is a blow-up of the boxed region in the left diagram. Black (colored) points denote saddle- 
node (grazing-sliding) bifurcations. Solid (dashed) curves represent stable (unstable) solutions. Lm 
decreases with increasing greenhouse gases. Horizontal shading around each transition in (b)-(c) 
indicates the amount of the smoothing (shading width is 2AE). 


diagrams using min(ii') instead of values chosen from a single Poincare map. 

3. We calculate the bifurcation diagram numerically for the albedo-smoothed system in (b) 
and (c) by integrating a range of initial conditions for one year to generate a Poincare map 
En+i = E{En), as in [3 [8]; each zero of the function G{E) = E — E{E) gives an initial 
value for a periodic state, and these initial conditions are used to compute min(Fl). For the 
bifurcation diagram for the Filippov system, we use the analytical solutions calculated in 
Section [3l For the seasonal state, the melting and freezing times tm and tf, defined implicitly 
by (13.19p . (13.201) . are found numerically. 

4. There are two gaps in the Filippov bifurcation diagram in (a); these gaps originate when 
solutions limit to the sliding intervals, as described in Section [3l (Recall that solutions that 
involve sliding are nonunique; we do not depict their corresponding solution branches.) The 
termination points of the solution branches in (a) occur where solutions limit to the sliding 
intervals; they correspond to grazing-sliding bifurcations. 

5. The seasonal state of the Filippov system always comes into existence through grazing-sliding 
as an unstable branch. This follows from the Floquet multiplier ()3.22p . which diverges at the 
existence boundaries associated with solutions limiting to the sliding intervals. 


14 












6 . In (a) and (b), notice that there are two hysteresis loops: a tiny one that involves the transition 
between the perennially ice-covered solution branch and the seasonally ice-free branch, and a 
large one that involves the perennially ice-free branch. The smaller hysteresis loop is shown in 
detail in the blown-up region in (a). This hysteresis loop is easily smoothed out by increasing 
the amount of smoothing in the model and is not present in (c). The absence of this bifurcation 
associated with the loss of summer sea ice was a major result of EW09. 

7. Most qualitative aspects of the bifurcation diagram of the albedo-smoothed system are pre¬ 
served in the Filippov limit; notice that this resemblance increases as AF decreases, which 
is illustrated in the difference between (b) and (c). (See Appendix for a more thorough 
comparison of bifurcations of the Filippov system and the albedo-smoothed system.) 


4.2 Impact of sliding intervals 

In this section we explore how the widths of the sliding intervals impact the size of the hysteresis 
loops in the corresponding albedo-smoothed system. We start with the observation that, roughly 
speaking, a wide sliding interval leads to a wide gap in the bifurcation diagram, whereas a narrow 
sliding interval creates a narrow gap. We find that for some parameter sets, narrow gaps may be 
associated with small hysteresis loops that are easily removed when smoothing is introduced. Figure 
0 ] provides two side-by-side examples demonstrating how the widths of the sliding intervals impact 
the sizes of the bifurcation diagram gaps. We use these examples as motivation for considering how 
the relationship between F+ and F- impacts the bifurcation structure. 

We wish to explore how changing the relative widths of the spring and autumn sliding intervals 
affects the bifurcation diagram. For this purpose we find it more convenient to write F± given by 
(j4.2h in the following standard form: 

F±{t) = F± + F± cos(27rr — ip±), (4-3) 

where 

1 Ajn i Act, 

-(1 ± Aa)5a - TaCos(27r())), (4.4) 

F_sin('0_) = —L(j sin(27r(^). 


_ F± = 
F±cos{^jJ±) = 
F+ sin(V’+) = 


from which it follows that 


F± = ^(1 ± Aq;)^^^ -I- 2(1 ± Aa)SaLa cos(27r(/)) -|- 
Aif; = 'ijj^ — ' 0 - 


= —sgn(())) arccos 


/ (l - Aa^) Si + 2SaLa cos(27r0) + Ll\ 

I F,f_ )' 


(4.5) 


where the factor —sgn((p) ensures that we have inverted the cosine function correctly; note that 
A-ip = 0 if (?!) = 0. 

The form of F± in (14.3p has the property that each parameter change causes a transparent 
transformation of the periodic function. It suggests a way to sweep through the parameter space 
of the model that is different from the parameter studies performed by Fisenman in [8] and is more 
systematic for our investigation of the impact of the sliding intervals. Note that the inverse of this 
mapping is not always unique: each parameter set given by (j4.2|] has a corresponding set in (14.3p . 
but parameter sets in (14.3p may map to parameters in (14.2p which are not permitted due to physical 
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^ L„ = 1.220 

Lm = 1.000 




(c) 


(d) 


Figure 4: Illustration of the relationship between the sliding interval widths and the bifurcation 
diagram gaps. (a),(c) show F+(r) (solid curve) and F_(r) (dashed curve) for the Filippov system 
()2.inp . (l4.2p : F± = 0 lines for different values are indicated by horizontal dashed lines. (b),(d) 
show corresponding bifurcation diagrams. Parameters of ()4.2I) are as in Figure [2l except in (c),(d) 
we set Sa = 1.58, La = 0.85, and (j) = —0.20. The sliding interval corresponding to each bifurcation 
diagram gap is indicated by a solid line in (a),(c) for the Lm values indicated by vertical lines in 
(b),(d). 


restrictions, or they may map to two sets of parameters. See Appendix |B] for a description of the 
inverse mapping we use in our investigation. 

Now we can explore how parameters in F±. change the bifurcation diagram. See Figure E] for 
several examples of how varying parameters of (j4.3p in turn, (a) A'i/’, (b) F+, and (c) T_, affect the 
bifurcation diagram. We make the following observations about Figure [5j 

1. The parameter F+ only affects the seasonal and perennially ice-free solution branches, while 
F- affects only the seasonal and perennially ice-covered states. 

2. The phase shift Aijj = — in contrast to F±, only impacts the seasonal solution branches 

of the bifurcation diagram and is perhaps the most interesting to vary. For example. Figure 
[5](a) shows that increasing A'lp, which controls the temporal offset between the peaks of F+ 
and F-, causes the seasonally ice-free branch to shift left to higher values. This causes 
the higher-energy gap in the bifurcation diagram to become smaller, while the lower-energy 
gap becomes larger. ^ ^ 

3. The widths of the sliding intervals, as each parameter (A^jJ, F^, F-) is varied, are shown 
in Figure [5](d)-(f). Here, for each parameter set we use the median value of in each 
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A-tp F+ F_ 


(d) 


(e) 


(f) 


Figure 5: (a)-(c) Bifurcation diagrams of the Filippov system (j2.inj) . (l4.2j) . with branches which 
change as an additional parameter of (j4.3p is varied overlaid in shades of gray, corresponding to 
black, dark gray, medium gray, and light gray branches as each parameter is increased (unchanged 
branches retain their original color). Diagrams correspond to (a) Aip = —0.21 (default parameters), 
0.00, 0.30, and 0.51; (b) F\ = 1.45, 2.00, 2.64 (default), and 4.00; and (c) f1 = 0.80, 
(default), 2.00, and 2.50. Parameters used to generate each diagram for the value of A'lp or F± 
shown are given in Table [2] in Appendix [Bj Note: the black and medium gray seasonally ice-free 
branches in (a) are identical to the seasonal branches in Figure[lKb),(d), respectively, (d)-(f) Curves 
showing the widths of the Si and S 2 sliding intervals (solid and dashed curves, respectively) as each 
parameter is varied. Markers highlight the sliding interval widths for the parameter values used in 
the corresponding hgure above. 


bifurcation gap to determine the width of the corresponding sliding interval. Notice that the 
two sliding intervals switch their relative widths for larger values of Atp, so that 1521 > |5i|, 
while 1 52 1 < |5i| for all values of F± considered. As illustrated in Figure O this behavior is 
related to the bifurcation diagram gaps changing their relative widths as increases. 

As Alp increases in Figure [51(a), the transition involving the perennially ice-covered solution 
branch, with increasing greenhouse gases, becomes large. A tiny hysteresis loop may form instead, 
which involves the transition between the perennially ice-free branch and the seasonally ice-free 
branch. Figure El shows the effect of reintroducing smoothing into the albedo function on the 
behavior of the bifurcation diagram in this extreme case, where = 0.51. Here, the smaller 

hysteresis loop that gets smoothed out is the transition from the perennially ice-free state to the 
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(a) 


(b) 


(c) 


Figure 6: Bifurcation diagrams for (I2.4l) - ()2.5p . where F{t,E) is given by ()4.1I1 . with (a) AE = 0, 
(b) AE = 0.02, and (c) AE = 0.08, for Aip = 0.51, obtained for Sa = 1.82, La = 1.19, cj) = —0.28, 
and other default parameters as in Figure [21 


seasonally ice-free state. In effect, increasing A'ljj has moved the small hysteresis loop from the low 
energy, thicker maximum ice-thickness end of the seasonal branch to the high energy, thinner max¬ 
imum ice-thickness end. There is now a sudden transition from the perennially ice-covered branch 
to the perennially ice-free branch. The smooth transition to the seasonal state with increasing 
greenhouse gases is completely lost, although there is now a smooth transition to a seasonal state 
from the ice-free state as greenhouse gases are decreased. 

The example presented in Figure [U] illustrates how the relative widths of the sliding intervals 
might impact the bifurcation diagram. It falls under Scenario 3 of Eisenman’s classification, which 
is concerned with possible abrupt transitions under increases in greenhouse gases, starting with the 
present-day perennial ice state [8]. However, varying a combination of the parameters used above 
may lead to more variation in possible bifurcation scenarios. To this end, we consider how the size of 
the jump from the end of the ice-covered branch to the seasonally ice-free branch, A(min(£')), varies 
while changing a pair of parameters. This proxy for the significance of an abrupt transition is small 
for the default parameters (A(min(E)) ~ 0.04; see Figure [3ja), zoomed box), but A(min(E)) may 
be increased by both decreasing F+ and increasing FT from their default values, the combination 
of which is similar to a decrease in Aa, while also increasing A'lp. 

Figure [7Da) shows a plot of A(min(E)) as a function of Aijj and Aa. A relatively large value 
of A(min(E)), approximately 0.3, occurs at {A'ip,Aa) ~ (0.27,0.3). The bifurcation diagram for 
these parameters is given in (b) for the Filippov model and (c) for the albedo-smoothed model. We 
make the following remarks about Figure [7l 

1. The horizontal dashed lines in (b) indicate where solutions jump off the ice-covered branch and 
hit the seasonally ice-free branch. The distance between the horizontal lines is A(min(E)). 

2. For a large region in (a), A(min(ii')) is similar to its value when calculated using the default 
parameters, where A(min(E)) ~ 0.04. However, for high values of Aip, A(min(E)) can 
become much larger, with the largest values occurring close to the upper black boundary. 

3. In the white region for Aa > 0.3 in (a), ice-covered solutions jump directly to the ice-free 
branch instead of the seasonal branch, as in Figure [Oja). The black boundary along the 
edge of this region is where the higher-energy saddle-node bifurcation on the seasonal branch 
coincides with the value of the ice-covered grazing-sliding bifurcation. The white region 
for lower Aa values (Aa < 0.3) represents where there are attracting sliding intervals; the 
black boundary here is where tc = (see Figured]). 
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A(min(i?)) 



Figure 7: (a) A(min(i?)) as a function of Aif; and Aa; inset: F±(t) for the default parameters 
(top) and for the parameters used in (b),(c) (bottom), (b) and (c) show bifurcation diagrams 
for {A'ij;,Aa) = (0.27,0.3), obtained for Sa = 1.68, La = 1.03, and (j) = —0.24, in (a) the Filippov 
model and (b) in the albedo-smoothed model with AE = 0.08 (cf. the default parameters: Sa = 1.5, 
La = 0.73, (p = 0.15, and Aa = 0.43). The distance between horizontal lines in (b) gives A(min(£')) 
for (AV’, Aa) = (0.27,0.3). 


4. The default value of Aa = 0.43 is determined in [8] using a^ = 0.68 and ami = 0.2. A 
decrease in Aa, given by (12.71) . to 0.3 may be achieved, for example, using Oj = 0.57 and 
keeping ami = 0.2. 

5. Notice that the jump to the seasonal branch in (b), A(min(Fl)) ~ 0.3, is smaller than the 
jump from the seasonal branch to the ice-free branch in the bifurcation diagram using the 
default parameter set, which is approximately 1.0 (see Figure [31(a)). 

6 . In (b), the jump to the seasonal state is much larger than the jump to the ice-free state. 
Note that there is no jump from the seasonal to the ice-free branch when the smoothing is 
introduced in (c). 

5 Discussion 

In our analysis we consider an energy balance model of the type introduced in EW09. The model 
in EW09 shows a smooth transition from a perennially ice-covered Arctic state to a seasonally 
ice-free state as greenhouse gases increase (see Eigure [31(c)). Through their analysis Eisenman and 
Wettlaufer showed that sea ice thermodynamics may play a key role in this smooth transition by 
counteracting the effects of positive ice-albedo feedback. In [8], Eisenman presented results of a 
thorough numerical exploration of bifurcation behavior in parameter space for a similar model. 


19 













determining the values where bifurcations and transitions between solution states occur as each 
parameter in (j4.ip is varied in turn, while fixing the others to their default values. He identified four 
possible bifurcation scenarios that the model produces, and discussed their physical interpretations 
and connections to results of other previously studied models. 

Our findings provide a case study that demonstrates how changing the relative widths of the 
sliding interval in a Filippov limit can affect hysteresis loops in bifurcation diagrams associated with 
the albedo-smoothed problem. Because we find that certain qualitative aspects of the bifurcations 
in the Filippov system are similar to their counterparts in the albedo-smoothed system, changes 
to the bifurcation diagram caused by varying the sliding interval widths are often reflected in the 
albedo-smoothed system as well. 

One uncertainty in these models, and in similar single-column models of the Arctic, lies in the 
representation of the albedo transition from that of ice to open ocean. This transition is often 
modeled as a sigmoid (see ( 12 .3p ). where the smooth transition parameterizes the spatial extent of 
the patchy region of ice between the ice-covered pole and open ocean at lower latitudes, as well 
as the dependence of albedo on bare ice thickness. However, model results are sensitive to the 
amount of smoothing. Additionally, a recent study by Wagner and Eisenman m investigated the 
limitations of parameterizations like the sigmoidal albedo transition in single-column models, by 
combining the EW09 single-column model with a spatially varying model. 

Our case study has a similar aim as previous studies of models of the type introduced in EW09: 
to understand how bifurcation behavior changes through parameter space. It is motivated by the 
relationship between the sliding interval widths and bifurcation diagram gaps, which are tied to the 
extent of hysteresis. We explore bifurcation behavior of the Filippov system by taking various paths 
through Eisenman’s parameter space which were suggested by the form of the periodic forcings and 
how they affected the widths of the sliding intervals. This provides an alternative perspective, 
where the focus is on the relative phase of the seasonal forcing (for an ice-free versus an ice-covered 
Arctic) in the presence of sea ice self-insulation feedback, rather than just on the role of the self¬ 
insulation feedback. From this perspective, we have identihed a scenario which may not be reflected 
in the scenarios considered in [ 8 ], e.g., the example in Figure [71(c). 

However, the relevancy of the example bifurcation diagrams given in Section 14.21 must be con¬ 
sidered. The parameters used in Figures [ 6 ] and [7] are within the range considered in the parameter 
studies performed in [ 8 ]. The largest difference in parameter values between the default set and 
those used in Figures [ 6 ] and [7] is the phase shift, 4>. The phase shift is the difference between when 
the maximum incoming shortwave radiation in ( 2 . 6 )) and minimum outgoing longwave radia¬ 
tion occur (J 7 in (2.6)). (Note that the minimum of J 7 , a loss term, corresponds to the maximum 
of the longwave forcing.) See Figure | 8 l(a) for an example of J-g^i using the seasonal forcing from 
the model in EW09 [3, nondimensionalized as in Eisenman’s analysis [ 8 ]. Figure (SKb) shows J-g^i 
from (14.2p using the default parameters, as in Figure [2] {4> = 0.15). Note the close correspondence 
between the forms of iFs,i in (a) and (b); this is a result of using a sinusoidal approximation of the 
observation-derived forcing iFs,i in (a) to determine the form and default parameters of J-g^i in (b) 
[ 8 ]. Figure ISKc) also shows from (|4.2I) . but using the parameter set corresponding to Aijj = 0.51 
instead, as in Figure El For (j) = —0.28, the phase shift in (c) (as with (j) = —0.24, the phase shift 
in Figure [7]), iFi reaches a minimum before the summer solstice, or alternatively, it is much delayed 
from the summer (due to the periodicity of the forcing). 

Physically, the seasonally varying component of the longwave forcing in Figure | 8 l(c) might 
occur in a different climate than the present-day (perennially ice-covered) state represented by the 
default parameter sets in both the model in EW09 [7] and the subsequent model by Eisenman 
[ 8 ]. There are two main contributions to the seasonal variation in the outgoing longwave radiation 
modeled in EW09: one contribution captures any seasonal variation in Arctic cloudiness, and the 
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Figure 8: Incoming shortwave radiation (solid curve) and outgoing longwave radiation 

(dashed curve), as defined by (12.611 with J) _|_ = Ti-- (a) From EW09 [7j, where J-g and J-i have 
been nondimensionalized by the annual mean incoming shortwave radiation and the mean absorbed 
shortwave radiation, respectively [8]. (b) From the sinusoidal approximation of EW09 used in [8], 
where Ts{t) = 1 — Sa cos ( 27 rT) and + La cos ( 27 r {t — cp)) from (j4.2p . using the default 

parameters, as in Eigure[2j (c) From (14.211 . using parameters, as in Figure^ that lead to Aip = 0.51. 


other contribution arises from the seasonal variation of atmospheric heat transport into the Arctic 
from lower latitudes [7]. In an alternative scenario, with perennially ice-free conditions, it has 
been suggested that a heavy cloud cover might form in the winter [2]. If this contribution were 
to dominate the seasonal variation of the outgoing longwave radiation, then it might shift the 
minimum of J-i from summer toward winter, as in Figure [81(c). 

Mathematically, this study provides an example of how bifurcation behavior of a smoothed 
system may be reflected in its Filippov counterpart, derived by taking the limit as a smoothing 
parameter goes to zero, and how the Filippov system may in turn help us understand behavior of 
the smoothed system. It also gives a unique perspective on the role of repelling sliding intervals in 
bifurcation behavior of a system with positive feedback that is periodically forced. The link between 
the width of the sliding intervals and the extent of hysteresis in the system is essential in exploring 
system behavior. We expect that this link may carry over to other systems with positive feedback, 
where the feedback strength is periodically modulated. Specifically, we anticipate the perspective 
used in this case study to apply when the positive feedback involves a sigmoidal function that can 
be taken to a PWS limit, as we did with the ice-albedo feedback. 


Appendix A Bifurcation set 

In this appendix we further explore how the bifurcation behavior of the Filippov system (I2.10|l . (l4.2p 
compares to that of the original albedo-smoothed system ()2.4l) . (j2.5l) . (l4.Ill used in the case study 
of Section |4l First we define the bifurcation points for the Filippov system using the existence 
conditions computed in Section [3l Then we show the bifurcation sets of both systems in the 
{Aa,Lm)- and ((/>, Lm)-parameter planes, holding all other parameters fixed at the default values, 
as given in Figure [2j We compare these results with corresponding figures given in Figure 9(F),(G) 
in [8]. We show that bifurcation behavior of the Filippov system is qualitatively similar to that of 
the albedo-smoothed system for the parameters studied. 

For the albedo-smoothed system investigated by Eisenman the only bifurcations identified in 
[8] are saddle-node bifurcations. However, he also keeps track of the existence regions for the 
three types of periodic states; for the albedo-smoothed problem, these existence boundaries are 
associated with the E = 0 boundary: perennially ice-covered (or ice-free) states become seasonally 
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Figure 9: Schematic bifurcation diagram illustrating the locations of the grazing-sliding points Lo,i 
and the saddle-node points Lsn^ 2 - 


ice-free once their maximum (or minimum) first hits FI = 0. In the context of nonsmooth systems, 
solutions which cross a discontinuity boundary as a parameter is changed are no longer “piecewise- 
topologically equivalent” to solutions that remain on one side of the boundary [6]. The parameter 
value where this collision with the FI = 0 discontinuity first occurs is a “grazing” bifurcation point 
of the albedo-smoothed system [6]. For the Filippov system, the existence boundaries are associated 
with grazing-sliding bifurcations, and there is also the possibility of saddle-node bifurcations on the 
seasonal solution branches. We label the grazing-sliding bifurcation point of the perennially ice-free 
(ocean) branch as Lq, the grazing-sliding bifurcation point of the perennially ice-covered branch as 
Lj, and the saddle-node bifurcation with the highest (lowest) minimum energy as (Tsn 2 ); see 
Figure m 

The values of Lq and Li are determined by the existence boundary for the corresponding peren¬ 
nial solution state, as derived in Section [3] f (l3.6j) and (I3.15h . respectively). The saddle-node bi¬ 
furcation points of the seasonally ice-free branch, simultaneously satisfy (13.231) and enable 

us to find which satisfy (j3.19|) . ()3.2nj) . Note that Lo^i always exist, but j col¬ 

lide and annihilate each other in a “hysteresis point.” The grazing-sliding bifurcations associated 
with the existence boundaries for the seasonally ice-free state have no direct counterpart in the 
albedo-smoothed system, so we omit them for our bifurcation curve study. Figure fTOl a) shows the 
bifurcation set for Fo,gsni 2 as a function of Aa, derived from our analysis in Section [3l The black 
region in (a) corresponds to where there is no repelling sliding interval, and the boundary for this 
region occurs where tj, = tc- 

Figure [H] presents a comparison of the bifurcation sets for the Filippov and albedo-smoothed 
systems in the (Aa, Lm)-parameter plane. In our analysis we use the grazing-sliding points Lo^i as 
proxies for the existence boundaries for the stable perennial branches, as suggested in Figure [3](b) 
and similar examples in Section 01 Figure [UK a) shows the bifurcation set for the values of Lo^i^sni 2 
as a function of Aa. In (b) we indicate the corresponding bifurcation set of the albedo-smoothed 
system, which is reproduced from [8] and was obtained in [8] by numerically constructing Poincare 
maps of the albedo-smoothed system. In (c) we overlay the bifurcation set of the Filippov system 
on the bifurcation set of the albedo-smoothed system for direct comparison. 

We make the following remarks about Figure [TTl 

1. The vertical line in (a) indicates the bifurcation points for the default parameters (cf. the 
bifurcation diagram of Figure [3K a)). 

2. The bifurcation set for the Filippov system in (a) shows the order and distance between 
bifurcation points as the additional parameter Aa is varied. 
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(b) 


(c) 


(d) 





(e) 


(f) 


(g) 


Figure 10: (a) Bifurcation set showing the values of Lo^i^sni ^2 (red, blue, and magenta curves, 
respectively) of the Filippov system (|2.1Up . (l4.2p as a function of Aa, using the default parameters 
in Figure [2j Vertical lines at Aa = 0.2, 0.3, and 0.43 show where bifurcation points occur in the 
bifurcation diagrams in (b),(c), and (d), respectively. (e),(f), and (g) show bifurcation diagrams of 
the albedo-smoothed system (|2.4l) . (j2.5l) . (|4.1l) for /S.E = 0.08, corresponding to the diagrams for the 
Filippov system directly above. 
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(b) 



0 0.2 0.4 0.6 0.8 1.0 
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Figure 11: (a) Bifurcation set for the Filippov system (|2.10p . (l4.2l) as a function of Aa, using the 
default parameters, (b) Bifurcation set and classes of stable fixed points of the albedo-smoothed 
system (iiaD,(i23|),(i3ii), reproduced with permission of the author from Figure 9(G) in [8]. Black 
boundaries between different classes represent smooth bifurcations; gray, red, and blue shading 
is where perennially ice-free, seasonally ice-free, and perennially ice-covered solutions are stable; 
light red is where both perennially and seasonally ice-free solutions are stable, and light blue is 
where both perennially ice-covered and ice-free solutions are stable (see [8] for further description), 
(c) Bifurcation set of the Filippov system overlaid on the bifurcation set of the albedo-smoothed 
system. 



Figure 12: (a) Bifurcation set showing the values of Lo^j^sni ,2 (red, blue, and magenta curves, respec¬ 
tively) for the Filippov system as a function of (j), using the default parameters, (b) Bifurcation set 
and classes of stable fixed points of the albedo-smoothed system, as in Figure 9(F) in [8]; shading 
is as described in Figure [TTl (c) Bifurcation set of the Filippov system overlaid on the bifurcation 
set of the albedo-smoothed system. 

3. The black region in (a) is similar to the white region of the albedo-smoothed system in (b), 
which for the smoothed system indicates where there is unphysical runaway cooling [8]. 

4. In (c), notice that the bifurcation set for the Filippov system is qualitatively similar to that 
of the albedo-smoothed system. In particular, note that Lo,i (the blue and red curves) are 
very close to the curves that signal the end of the perennial solution branches in the albedo- 
smoothed system, as illustrated in Figure [3Kb). 

In Figure [121 we provide a similar comparison of the bifurcation set between the Filippov and 
albedo-smoothed systems, for varying </>. Note that our analysis requires us to exclude a region of 
parameter space not excluded in the albedo-smoothed system, due to the presence of attracting 
sliding intervals in the Filippov system. Also note that in (b) we have shifted the range of the 
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phase (p from the original figure in [ 8 ] so that p E (—0.5,0.5]. In general, while there are noticeable 
quantitative differences, there is significant qualitative similarity, as in Figure fTTl 


Appendix B Inverse mapping 

In this appendix we determine an inverse mapping, when it exists, between Eisenman’s physical 
parameters, Aa, Sa, Lm, La, p of (|4.2I) . and the parameters we introduced for our analysis in 
Section I121 namely, F±, F± and Ap = p-^- — p- of (I4.3p . We are interested in the mapping so that, 
for example, we can construct a path through Eisenman’s five-parameter space (L^, Aa, Sa, La, 
p) that is associated with varying the relative phase Ap, while holding the four amplitudes (F+, 
F_, E+, E_) fixed. At the end of this appendix is a table of parameter sets used to generate the 
bifurcation diagrams in Figure [5l 

Here we reproduce the two expressions that are relevant to our construction: 

L'±{t) = (1 ± Aa)(l - Sa cos(27rr)) - (L^ + La cos(27r(r - p))) (B.l) 


and ^ 

F±{t) = F± + F± cos(27rr — p±). (B.2) 

We can relate the phases {'p±) and amplitudes {F±,F±) to Eisenman’s physical parameters via 


_ F± 
F± cos{p±) 
F±sm{p±) 


I Lrfa zb Acx, 

-(1 zb Aa)Sa - La cos(27r(/)), 
—La sin(27r(/)). 


(B.3) 


We introduce the following shorthand for some of the relationships in (IB.31) : 


where 


F+cos{p+) 
F_ cos(V’-) 
F±sm{p±) 


* 5 '+ Lac, 
^<5*+ Lac, 
Las, 


S+ 

r 

Lac 

Las 


{1 + Aa)Sa € [Sa,2Sa], 


1 — Aa 
1 -b Aa 


€ [ 0 , 1 ], 


LaCOSplTTp) E [-La, La], 
LaSinipiirp), and where L^^ + L\ 


Li 


Here we enumerate some preliminary steps used in constructing an inverse mapping: 
1. F+ and F- determine both Lm and Aa, which in turn determine r in (IB.4j) - ()B.5j) : 


A„ = i(F,-F_), 

Lm ~ 2 ('^+ ■ 


(B.4) 


(B.5) 


(B. 6 ) 


Note that the conditions Aa E [0,1] and L^ > 0 place some restrictions on the possible 
values of F±. 
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2. We next use the values of -F+ and F_ from (|B.4h to determine the parameter combinations 
5+ and Lac, where 

F+ = Ll + Sl + 2S+Lac, 

F? = Ll + r^Sl + 2rS+Lac. 

This pair of equations can be solved, given F± and r determined in step 1, to determine 5+ 
and Lac as functions of Specifically, we find 


^1 = 


S-i-Ln.r. — 


rF+^ - fJ + Ll{l - r) 


r(l — r) 

~2 ~2 

— F + La{l — 


2r{r — 1) 


(B.7) 


The positive root of the equation is used to determine 5+ > 0. Note that there are some 
restrictions placed on the possible values for F± by the physically allowed ranges on and 
Lac- 

3. We next obtain an equation for La in terms of the quantities calculated in the previous steps 
and tan(A'i/^), where 


tan(A'!/;) 


sin(?/;_i_) cos{'tp-) — cos{^p+) sin('0-) 
cos('0+) cos{'tp-) + sin('!/)+) sin(V’-) 
(r - 1)5+La^ 

La + r 5"^ + Sj^Laci^ + r) 


(B.8) 


The values of 5^ and S+Lac are given by ()B.7jl in terms of the amplitudes F± and 
Also, Las = where the choice of the sign will ultimately be determined by 

the sign of tan(AV’), and in turn that will determine the sign of cj) since Las = LaSm{2Tr4>), 
where (j) G Combining all of this, we find that solves the following quadratic in 

F+, F_, r, and tan^(A^/;): 


4 2iFJ + PF+] 

“ (l-r)2 


+ 


(T_ — )^ + (F_ + )^tan^(AV') 

(1 — r)^(l + tan^(AV^)) 


= 0 . 


(B.9) 


Again, there are restrictions on the values of F±. and A'lp in order for this to have real solutions. 
Moreover, there can be two positive, real solutions to this equation. We find, when there 
are two distinct solutions for that one is much larger than the other, and it is typically 
the smaller one that is of interest to us. 

To summarize, given a set of parameters F±, F± and A^jJ of (jB.2p . we compute parameters for 
the Eisenman model (IB.II) as follows: 

1. Compute Aa and Lm using (IB.611 . and then compute r using its definition in ()B.5jl . 

2. Compute La from the quadratic (IB.9p . which may have two roots. The choice of root deter¬ 
mines La, which is nonnegative. 

3. Compute 5+ and Lac from (IB.71) . 5+ and Aa determine Sa using (IB.51) . 

4. Las can be determined, up to a sign, from and Lac- Then the sign of Las is determined 
from the sign of tan(A?/:) in (IB.Sp . Finally the phase 4> is determined from Lac and Las, 
defined in (IB.Sp . 
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Mapped parameters 

Original parameters 


Alp 

F+ 

f1 

Sa 

La 



-0.21 

2.64 

1.41 

1.50 

0.73 

0.15 

Varying Aif 

0.00 

2.64 

1.41 

1.43 

0.60 

0.00 

0.30 

2.64 

1.41 

1.58 

0.85 

-0.20 


0.51 

2.64 

1.41 

1.82 

1.19 

-0.28 


-0.21 

1.45 

1.41 

0.34 

1.42 

0.27 

Varying F+ 

-0.21 

2.00 

1.41 

0.79 

1.09 

0.14 

-0.21 

2.64 

1.41 

1.50 

0.73 

0.15 


-0.21 

4.00 

1.41 

3.06 

0.59 

0.38 


-0.21 

2.64 

0.80 

2.17 

0.52 

0.43 

Varying E_ 

-0.21 

2.64 

1.41 

1.50 

0.73 

0.15 

-0.21 

2.64 

2.00 

0.93 

1.65 

0.15 


-0.21 

2.64 

2.50 

0.63 

2.47 

0.25 


Table 2: Approximate parameters from (14.2p used to generate the bifurcation diagrams in Figure O 
calculated using the inverse mapping described in this appendix. The remaining parameters of the 
full model {B, and Aa) have been set to the values given in Figure [2j 
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